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Abstract. The phylogeographic structure of some species distributed across the Baja California Peninsula has been 
traditionally hypothesized as resulting from vicariant events thought to have occurred between 1-3 Mya. Climatic 
fluctuations during the Pleistocene have also been shown to influence the distribution patterns of species, and vicariant 
patterns may have been erased as a consequence of population contractions or expansions into or out of refugia generated 
during the last glacial maximum ca. 21,000 years ago. Thus, there is still some uncertainty regarding the relative role of 
vicariance in shaping the modem biota of Baja California. To understand the evolutionary history of the wolf spider 
Pardosa sierra Banks, 1898 on the peninsula, a phylogeny of this species and closely related taxa was generated using a 
fragment of the mitochondrial gene cytochrome c oxidase subunit I (COl). Sequences of a fragment of the COl gene for 38 
individuals from 14 sampling sites along the entire distribution range of P. sierra were used to infer phylogeographic 
patterns, and five nuclear microsatellite loci were also used to genotype 296 individuals from seven of these 14 locations. 

The current and past potential distributions from two Pleistocene periods were estimated using niche-based distribution 
modeling, and scenarios of colonization from detected refugia were simulated. We found that Californian populations of 
P. sierra diverged from peninsular populations 4 Mya, this divergence coinciding with the northern-gulf split. However, we 
did not detect genetic breaks in regions where the mid-peninsular and Isthmus of La Paz canals were presumably formed, 
either with mitochondrial DNA sequences or microsatellite loci. Two refugia were further detected at the geographic ends 
of the peninsula, these likely preceding subsequent habitat expansion. 
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The Baja California Peninsula (BCP; Fig. 1) has been a 
model landscape for phylogeographic studies due to its 
geographic isolation, landscape heterogeneity and high levels 
of endemism (Jezkova et al. 2009; Wilson & Pitts 2012; 
Graham et al. 2014; Dolby et al. 2015). In several studies, 
population divergences or phylogeographic breaks have been 
detected in a diversity of taxa (Upton & Murphy 1997; Riddle 
et al. 2000; Murphy & Aguirre-Leon 2002; Nason et al. 2002; 
Zink 2002a; Crews & Hedin 2006; Lindell et al. 2006; Graham 
et al. 2014; Lira-Noriega et al. 2015). To explain the causes of 
such phylogeographic breaks, hypotheses that include vicar¬ 
iance events have been proposed, with the most important 
being these: (1) the marine transgression of the Gulf of 
California in the northern-gulf region (south of California and 
Arizona), with separation from the rest of the BCP occurring 3 
million years ago (Mya) (Riddle et al. 2000; Hafner & Riddle 
2005); (2) the mid-peninsular channel, which formed approx¬ 
imately 1 to 1.6 Mya in the nearby Vizcaino Desert (Upton & 
Murphy 1997; Hafner & Riddle 2005; Crews & Hedin 2006; 
Lindell et al. 2006); and (3) the Isthmus of La Paz Channel, 
which separated the Cape Region of the BCP 3 Mya (Riddle et 
al. 2000; Hafner & Riddle 2005; Garrick et al. 2013). 

Alternative non-vicariant biogeographic hypotheses have 
also been proposed, most of these highlighting the potential 
effects of cyclic climatic changes during Pleistocene glacial and 


interglacial periods (Hewitt 1996, 2004; Hafner & Riddle 2005; 
Lindell et al. 2006; Riddle & Hafner 2006; Leache & Mulcahy 
2007; Garrick et al. 2009, 2013). In particular, the Last Glacial 
Maximum (LGM) ca. 21,000 years ago (Dansgaard et al. 
1993) may have forced many species from the deserts of the 
Northern Hemisphere to find refuge south of their previous 
distributions (Hewitt 1996, 2000, 2004; Van Devender 2002; 
Graham et al. 2013; Harl et al. 2014). The Cape Region, at the 
southern tip of the BCP, was likely used as a refugium for 
species that could not tolerate cold (Murphy & Aguirre-Leon 
2002; Garrick et al. 2013). However, the number and location 
of refugia is species dependent, as each likely responded 
differently to climatic events (Harl et al. 2014). Similarly, it is 
also probable that cycles of range contraction and range 
expansion generated by climatic oscillations (Grismer 2002) 
may have erased the genetic fingerprints produced by earlier 
vicariance events (Crews & Hedin 2006). This raises uncer¬ 
tainty about whether the phylogeographic patterns observed 
are due to hypothesized vicariance events or climatic changes 
or both. 

Pardosa sierra Banks, 1898 is a species of wolf spider 
(Lycosidae) that is endemic to the BCP (Correa-Ramirez et al. 
2010); it is intolerant of cold temperatures (Van Dyke & 
Lowrie 1975) and strongly dependent on water bodies (Punzo 
& Farmer 2006; Correa-Ramirez 2010; Jimenez et al. 2015). 
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Figure 1 .—Pardosa sierra collection sites in the Baja California 
Peninsula. For the acronyms of each collecting site refer to Table 1. 
Lines indicate the hypothetical vicariance events of the northern-gulf 
{blue), the Vizcaino mid-peninsular channel (red) and the Isthmus of 
La Paz (green). CA = California; BC = Baja California; BCS = South 
Baja California; LP = La Paz. 

Given these biological constraints, it is expected that 
populations of P. sierra would have been affected by changes 
to their habitat due to climatic variations during the 
Pleistocene, as has been shown for other desert species of 
Rodentia (Jezkova et al. 2009), Hymenoptera (Wilson & Pitts 
2012), and Scorpiones (Graham et al. 2013), among others 
(Van Devender 2002). 

In this study, mitochondrial DNA sequences and nuclear 
microsatellites were used to analyze the phylogeographic 
patterns and genetic structure of P. sierra populations across 
the BCP, with the aim of addressing the following questions: 
What is the phylogeographic pattern of P. sierra along the 
BCP? Flow and when did current distribution patterns 
potentially originate? Does the genetic structure inferred using 
microsatellite loci coincide with phylogeographic patterns 
detected using mtDNA sequences? Did P. sierra use refugia in 
more than one climatically suitable area during the LGM? If 
so, where were these refugia? We expect that any genetic signal 
resulting from potential vicariance events that occurred before 
or during the Pleistocene were erased due to refuge- 
colonization processes caused by climatic changes after the 


LGM. Finally, if P. sierra expanded its distribution from 
southern refugia after the LGM, we predict a decreasing 
gradient of microsatellite genetic diversity from the south 
region to the north of the BCP. 

METHODS 

Samples and study locations. —Between 2006 and 2011, 296 
adult P. sierra specimens were hand collected from the 
margins of rivers, ponds and other suitable habitats along 
the BCP and southwestern United States (Fig. 1; Table 1). We 
collected 38 specimens from 14 localities throughout the range 
of the species for mitochondrial DNA sequence analysis, and 
296 specimens were sampled from only seven locations along 
the BCP for microsatellite genotyping (Table 1). The 
specimens were preserved in 95% ethanol and deposited in 
the Coleccion Aracnologiea y Entomologica del Centro de 
Investigaciones Biologicas del Noroeste (CAECIB). 

DNA extraction and genetic characterization using mito¬ 
chondrial markers and microsatellites. —Total genomic DNA 
was extracted from the legs of each specimen as reported by 
Aljanabi & Martinez (1997). DNA quality was verified on 1% 
agarose gels (buffer TAE lx. Gel Red lOx), and DNA 
concentration (ng/pl) was quantified with a NanoDrop 8000 
(Thermo Fisher Scientific, Wilmington, DE). A mitochondiral 
fragment of the cytochrome c oxidase subunit I (COl) gene 
was amplified by polymerase chain reaction (PCR) and 
sequenced using the primers CO IP- L and COIP-R (Correa- 
Ramirez et al. 2010). PCR reactions were performed with a 
PTC-200 DNA Engine Thermal Cycler (BioRad Laboratories, 
Hercules, CA) in a total reaction volume of 15 pi (~50 ng 
genomic DNA, 0.40 mM each primer, 2.5 mM mM MgCl 2 , 0.2 
mM each dNTP, lx PCR buffer and 0.5 U Taq polymerase 
(Invitrogen, Carlsbad, CA)). A total of 35 temperature cycles 
were performed, each cycle consisted of denaturation for 30 
seconds at 94°C, annealing for 30 seconds at 52°C, extension 
for 30 seconds at 72°C, and included an initial denaturation 
step for 4 minutes at 94°C, and a final extension step for 10 
minutes at 72°C. PCR products were visualized using 
electrophoresis in 1.5% agarose gels. The sequencing of PCR 
products was performed using the BigDye Terminator 
sequencing method in an ABI PRISM 3730XL sequencer 
(PE Applied Biosystems; Macrogen, Seoul, Korea). 

Five microsatellite loci previously described for P. sierra in 
Molecular Ecology Resources Primer Development Consor¬ 
tium et al. (2010) were used. The final PCR reaction volume 
was 15 pi with ~30 ng genomic DNA, lx PCR buffer (20 mM 
Tris-HCl, pH 8.4, 50 mM KC1), 200 pM each dNTP, 0.4 pM 
each primer, 1.5-2.5 mM MgC12 and 0.15 U Taq DNA 
Polymerase (Invitrogen, Carlsbad, CA). The PCR conditions 
were 95°C for 5 min, followed by 35 cycles of 95°C for 30 see, 
annealing at 52°C for 30 sec, 72°C for 30 sec, and a final 
extension of 72°C for 5 min. PCR reactions were performed in 
an MJ Research PTC-200 thermal cycler. PCR products were 
visulaized by electrophoresis on 6% polyacrylamide gels with 
7.5 M urea. Gels were silver-stained according to Benbouza et 
al. (2006). The allele size was determined by their relative 
position on the gel compared to the molecular marker 
(nucleotide ladders of 10 and 50 bp, Invitrogen). 

Phylogeography and population genetic structure,, - To assess 
whether the phylogeographic pattern of P. sierra along the 
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Table 1.—Sampling locations of Pardosa sierra in the Baja California Peninsula and frequency of haplotypes of mitochondrial gene COl. if = 
number of individuals analyzed for each microsatellite locus. 


Location 

ID 

Lat (°N) 

Long (°W) 

/? a 

Haplotypes COl (n) 

USA, CA, Sn. Juan Creek 

SNJC 

33.638 

-117.421 

_ 

Hap3 (4) 

USA, CA, Sn. Felipe 

SNFE 

33.066 

-116.553 

- 

Hapl (1), Hapl (3) 

Mexico, BC, Arroyo Las Palomas 

ARLP 

32.373 

-116.354 

- 

Hap\ (2), Hap6 (2) 

Mexico, BC, Arroyo Sn. Antonio Minas 

ASNM 

31.968 

-116.658 

- 

Hap 1 (4) 

Mexico, BC, Ensenada 

ENSE 

31.783 

-116.6 

26 

Hap 1 (1), Hap5 (3) 

Mexico, BC, Rancho Las Liebres 

RLL1 

31.584 

-116.03 

- 

Hapl (2) 

Mexico, BC, Arroyo El Mejin 

AREM 

30.980 

-116.094 

- 

Hap6 (1) 

Mexico, BC, El Rosarito 

ROSA 

28.616 

-114.033 

10 

Hap 1 (1), Hap5 (1) 

Mexico, BCS, Sn. Ignacio 

SNIG 

27.174 

-112.869 

- 

Hap 1 (1) 

Mexico, BCS, Cadeje 

CADE 

26.366 

-112.5 

52 

Hapl (2), Hap4 (2) 

Mexico, BCS, Sn. Isidro-La Purfsima 

SNIP 

26.2 

-112.033 

52 

Hapl (1) 

Mexico, BCS, Sn. Pedro de la Presa 

SNPP 

24.833 

-110.983 

52 

- 

Mexico, BCS, Presa de la Buena Mujer 

PBMU 

24.088 

-110.191 

- 

Hapl (2) 

Mexico, BCS, El Novillo 

NOVI 

23.916 

-110.216 

52 

Hapl (3) 

Mexico, BCS, Sierra de la Laguna 

SLLA 

23.233 

-109.866 

52 

Hapl (2) 


Baja California peninsula was correlated with geological 
vicariance events, the 14 locations were clustered a priori into 
four groups, California (CA), Baja California (BC), Baja 
California Sur (BCS) and La Paz (LP), based on the positions 
of the hypothetical trans-peninsular channels (Fig. 1). An 
analysis of molecular variance (AMOVA) was performed 
using COl sequences among these groups. Additionally, 
genetic differentiation values ( F S t ) were estimated by applying 
10,000 permutations using ARLEQUIN 3.5 (Excoffier & 
Lischer 2010). 

To determine the genetic structure of P. sierra . microsatel¬ 
lite loci were analyzed using STRUCTURE 2.3.4 software 
(Pritchard et al. 2000). First, we performed an analysis without 
the LOCPRIOR option to obtain an initial K , and afterwards 
we ran it again using this K value as LOCPRIOR. In both 
analyses the admixture model was used, with a burn-in period 
of 100,000 and lxlO 6 MCMC repetitions thereafter. The 
number of populations was estimated using the A K value using 
Evanno’s method (Evanno et al. 2005) in STRUCTURE 
HARVESTER Web vO.6.94 (Earl & VonHoldt, 2012). The 
results of 15 independent runs were processed and visualized 
using CLUMPP 1.1.2 (Jakobsson & Rosenberg 2007) and 
DISTRUCT 1.1 (Rosenberg 2004) respectively. Finally, the 
genetic differentiation values ( F S t ) were evaluated using the 
infinite allele model (Weir & Cockerham 1984) in ARLE¬ 
QUIN, using the Bonferroni correction to adjust the 
significance level for multiple tests (Rice 1989). 

Genetic diversity.—To observe haplotype diversity, haplo- 
type frequencies from the COl sequences were calculated 
using DnaSp 5.10 (Librado & Rozas 2009). To estimate 
population growth of P. sierra , the R2 test was performed 
(Ramos-Onsins & Rozas 2002), with 10,000 coalescence 
simulations using DnaSP v5 (Librado & Rozas 2009). To 
detect potential relationships among haplotypes, a statistical 
parsimony network was estimated for COl with a 95% 
confidence criterion, using TCS 1.13 software (Clement et al. 
2000 ). 

To detect a genetic diversity gradient along the peninsula, 
we calculated the Pearson correlation coefficient between the 
latitude of each location with its corresponding heterozygosity 
(. H e ) and allelic diversity (A/L) of microsatellites previously 


obtained from GENEPOP 4.0.7. (Rousset 2008). To observe if 
there is a relationship between distance and genetic identity 
between populations, we tested for isolation-by-distance, using 
the Mantel test (Slatkin 1993) with 10,000 permutations in 
IBDWSin v3.23 (Jensen et al. 2005). Recent population 
bottleneck analysis was performed with a two-tailed Wilcoxon 
sign-rank test for heterozygosity excess under a two-phased 
mutation model (TPM; Di Rienzo et al. 1994; Miller et al. 
2012), using the program BOTTLENECK 1.2.02 (Piry et al. 
1999). 

Phylogenetic analysis.—In addition to the newly-generated 
COl sequences from this study, additional sequences were 
obtained from GenBank (http://www.ncbi.nlm.nih.gov/ 
genbank/; accessed April 2014) for another 28 taxa, these 
included different species of Pardosa C.L. Koch, 1847 and 
closely related genera (Appendix 1). For Bayesian phyloge¬ 
netic inference we used MR BAYES 3.2.2 (Ronquist & 
Huelsenbeck 2003), and JMODELTEST 2.1.6 software 
(Darriba et al. 2012) was used to identify an optimal 
GTR+I-bG model of molecular evolution for the un-parti- 
tioned COl matrix under the Akaike Information Criterion 
(AIC). Bayesian analysis was performed using two indepen¬ 
dent runs consisting of 10xl0 6 MCMC generations, with 
sampling every 1,000 generations. A 50% majority-rule 
consensus tree was generated after 10% burn-in. 

Estimation of divergence time.—To test whether phylogeo- 
graphic patterns detected were temporally concordant with 
hypothetical vicariance events, interspecific and population 
divergence times were estimated in BEAST 1.8 (Drummond et 
al. 2012), using both fossil and rate calibration methods with 
COl. The first analysis included individuals of P. sierra , plus 
other Pardosa species and related genera. The GTR+I4G 
substitution model was used based on the results of 
JMODELTEST 2.1.6, along with a lognormal relaxed clock 
model and Yule speciation process to model the tree prior. 
The Yule process evaluates the relative ages of nodes 
contributing to the prior distribution of nodal ages (Ho & 
Phillips 2009). A lycosid macrofossil from the Miocene found 
in Dominican amber (Penney 2001) was used as a root 
calibration point, with a minimum age constraint of 20 Mya. 
We implemented a hard minimum bound (Ho & Phillips 2009) 
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to a well-identified macrofossil with a lognormal distribution 
(mean 0, SD 1.0, offset 20). The fossil specimen used here is 
considered the first representative of the family Lycosidae 
from the fossil record (Penney 2001); previously recorded 
fossils from the Carboniferous, Tertiary, Miocene, or Pliocene 
assigned to genera of Lycosidae were likely misidentified 
species that correspond to other families (Penney 2001). 

We also performed a second BEAST analysis including only 
P. sierra specimens, using the Coalescent constant size to 
model the tree prior. The tree was rate calibrated using the 
substitution rate of 0.0115 substitutions per site per million 
years for COl. proposed for insects by Brower (1994) and also 
used for spiders in phylogeographic studies (e.g., Chang et al. 
2007). For both BEAST analyses, two independent runs of 
40x10 6 generations were performed, sampling every 2,000 
generations. The results were analyzed in TRACER 1.5 
(Rambaut & Drummond 2007). To ensure parameter conver¬ 
gence and the effective sampling size (ESS), all parameters and 
trees for both independent runs were combined in LOG- 
COMBINER 1.8 (Rambaut & Drummond 2007) with a burn- 
in of 25% for the first trees in TREEANNOTATOR 1.8 
(Rambaut & Drummond 2007). Finally, the results were 
summed up in a single tree (maximum clade credibility tree) 
and visualized in FIGTREE 1.4 (http://tree.bio.ed.ac.uk/ 
software/figtree/). 

Niche-based distribution model. To identify geographic 
areas that might have served as refugia for P. sierra during 
and after the LGM, we generated a niche-based distribution 
model (NBDM; Segurado et al. 2006) using records of P. 
sierra collected for this study at 31 localities in the States of 
Baja California and Baja California Sur in Mexico, and 
California in the United States of America (Appendix 2). To 
infer the NBDMs the records were loaded into the 'maximum 
entropy machine-learning algorithm'’ MAXENT 3.2.2 (Phil¬ 
lips et al. 2006). Current bioclimatic variables used (BIO 1-19) 
were downloaded from WorldClim (Hijmans et al. 2005; www. 
worldclim.org) with a resolution of 1 km 2 using 75% of the 
records for training and the remaining 25% for validation of 
the model. To run MAXENT, default parameters were used 
with 1,000 iterations. The models were evaluated with the area 
under the curve (AUC) method. The algorithm compensates 
for co-linearity between variables using a method for 
regularization that deals with feature selection. There is thus 
less of a need to remove correlated variables (Elith et al. 2011), 
as the algorithm ranks the contribution of each during the 
analysis. The results were projected in QUANTUM GIS 2.2 
Valmiera software. To explore the potential occurrence of a 
species in the past, the models generated under current 
climatic conditions were projected onto the LGM (21,000 
years ago) and interglacial (140,000 -120,000 years ago) 
scenarios. The climatic layers from the past were downloaded 
from WorldClim for LGM scenarios developed by the 
Paleoclimate Modelling Intercomparison Project Phase II 
(Braconnot et al. 2007), the Community Climate System 
Model (CCSM: Collins et al. 2004), the Model for Interdis¬ 
ciplinary Research on Climate (MIROC; Hasumi & Emori 
2004) and, for the interglacial period, we used the layers 
developed by Otto-Bliesner et al. (2006). The CCSM and 
MI ROC models simulate climatic conditions that were 
calculated as prevailing during the LGM, however a stronger 


decrease in temperature is assumed for the CCSM compared 
to the MI ROC model (Otto-Bliesner et al. 2006). For assessing 
variable importance in each model a jackknife test was 
performed in MAXENT. 

RESULTS 

Phylogeography and population genetic structure.—The 

mitochondrial COl gene studied had 630 base pairs, and in 
the sample of 38 specimens, 7 unique haplotypes were detected 
(Fig. 2a, b; Table 1; GenBank accession numbers in Appendix 
1). No population expansion was detected with the COl gene 
(R2 = 0.099; P= 0.239). The haplotype network shows that the 
haplotypes tend to be grouped according to their geographical 
regions. Haplotype 1 (Hap]) was the most abundant in the 
southern region of the BCP, whereas Hap2 and Hap6 were 
only detected in the northern region of the BCP. Additionally, 
Hap3 and Hap? were only observed in California (USA), 
seven mutational steps apart with respect to the remaining 
haplotypes (Fig. 2a, b; Table 1); this coincides with the 
geographical separation of clades detected in the phylogenetic 
analysis (Fig. 3). The AMOVA detected that the Californian 
locations are different from all the BCP locations (F cr =Q.7, P 
= 0.001), whereas inside the BCP using another test, no 
differences were observed between the BC, BCS and LP 
groups ( F C t= 0.1, P — 0.07) (results not shown). STRUC¬ 
TURE clustered the seven locations of the BCP into five 
populations, although genetic mixing exists between them 
(Fig. 2 c, d), which in north-south latitudinal order are I 
(Ensenada), II (El Rosarito and Cadeje), III (San Isidro La 
Purisima and El Novillo), IV (San Pedro de La Presa) and V 
(Sierra de La Laguna). Pairwise F$ T values showed genetic 
differences between BCP populations (F sr varied from 0.01 to 

O. 04, P < 0.001, Table 2). 

Using microsatellite data, the genetic diversity detected was 
9.4-17.6 for A/L and between 0.81-0.84 for II E (Table 2). The 
Mantel test for the microsatellites (/• = 0.36, P = 0.9), and 
correlation analyses between latitude and the population 
heterozygosity (/• = -0.64, P = 0.2), and between latitude and 
A/L (r = —0.82, P = 0.08) were not significant. The bottleneck 
analysis was significant for the Wilcoxon sign-rank test only 
for northern populations I (P = 0.03) and II (P = 0.03). 

Phylogenetic relationships.—The Bayesian 50% majority- 
rule consensus tree (Fig. 3) recovered the samples of P. sierra 
as monophyletic (PP= 1.0), and sister to P. atromedia Banks, 
1904 (but with low nodal support; PP = 0.69). The sister- 
groups to P. sierra + P. atromedia were a polytomic 
assemblage of taxa including P. valens Barnes, 1959, P. steva 
Lowrie & Gertsch, 1955, P. sura Chamberlin & I vie, 1941, and 

P. vadosa Barnes, 1959. Within the P. sierra lineage, only two 
subclades were detected: the California clade (CA) from the 
USA, and the Baja California Peninsula (BCP) clade from 
Mexico, each reciprocally monophyletic. This analysis did not 
detect divergence among populations of P. sierra from the 
BCP (Fig. 3). 

Estimation of divergence times.—BEAST analysis using a 
fossil calibration inferred a divergence date for the separation 
of the CA clade from the BCP clade of ca. 4 Mya (95% HPD 
8.52-1.65 Mya; Fig. 4), The CA clade started diversifying 
around 1.6 Mya (95% HDP 5.9-0.76 Mya) and the BCP clade 
around 2.4 Mya (95% HDP 4.3-0.3 Ma). Alternatively, for the 



Figure 2. -Fhylogeographic pattern and genetic structure of Pardosa sierra in the Baja California Peninsula, (a) Statistical parsimony network based on 38 COl sequences. The areas 
of the box and circles are proportional to the haplotype frequency, and each nucleotide substitution between haplotypes is represented by a dot. (b) Latitudinal distribution of seven 
haplotypes of the mtDNA gene COl. The acronyms of each collecting site are described in Table 1. (c) Number of genetically homogeneous populations (I V) detected, using five 
microsatellite loci in STRUCTURE. Abbreviation (n) indicates sample size per population, (d) In the diagram, each horizontal bar represents an individual and each color represents the 
membership proportion (qi) corresponding to each population defined by STRUCTURE. Lines and arrows on the maps indicate the putative channel transgressions of the northern-gulf 
(blue in b), the mid-peninsular Vizcaino channel (red in b, c, and d) and the Isthmus of La Paz (green in b, c, and d). 
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Figure 3.—Bayesian 50% majority-rule consensus tree of mitochondrial COl sequences, for populations of Pardosa sierra and related 
Lycosidae. Numbers below the nodes indicate the posterior probability values. BCP = Baja California Peninsula, MEX = Mexico, CA = 
California, USA. For acronyms of the location of each specimen see Table 1. 


Table 2.—Values of genetic differentiation ( F ST ) estimated by pairs 
of populations, allelic diversity (A/L), expected (H E ) and observed 
(Ho) heterozygosity, based on the allelic frequencies of five 
microsatellite loci of Pardosa sierra from the Baja California 
Peninsula. Values below the diagonal show the pairwise F S r values 
and above the diagonal are shown the corresponding probability 
(values *** P < 0.0001). I (Ensenada), II (Rosarito-Cadeje), III (San 
Isidro La Purisima and El Novillo), IV (San Pedro de la Presa), V 
(Sierra de La Laguna). 




Locations 



A/L 

Ho 

He 

I 

II 

III 

IV 

V 

I 

_ 

*** 

*** 

*** 

*** 

9.4 

0.55 

0.81 

II 
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rate-calibrated BEAST analysis, the split between both clades 
(CA and BCP) was inferred as being much more recent (0.6 
Mya; 95% HPD 1.04-0.22 Mya; tree not shown), with 
diversification from 0.1 Mya (95% HDP 0.23-0.007 Ma) and 
0.2 Mya (95% HOP 0.37-0.06 Ma) for the CA and BCP 
clades, respectively. 

Niche-based distribution modeling.—The model projected for 
the interglacial period (140,000-120,000 years ago) indicated 
expanded areas of potential habitat in the northern and 
southern BCP (Fig. 5a). The models CCSM and MI ROC 
projected for the LGM (21,000 years ago) revealed that the 
potential habitat of P. sierra was concentrated in the 
latitudinal geographic ends of the BCP and that an important 
wide band of inadequate habitat was projected in the center of 
the BCP (AUC = 0.75-0.99 for both models; Fig. 5b. c). 
Currently, the distribution predicted by MAXENT agrees 
with the known distribution of P. sierra (AUC = 0.945; Fig. 
5d). The most important climatic variables identified by the 
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Figure 4.—Fossil calibrated BEAST chronogram of mitochondrial COl sequences, for populations of Pardosa sierra and related Lycosidae. 
Green bars indicate the 95% highest posterior density estimates of divergence time for each node. BCP = Baja California Peninsula, MEX = 
Mexico, CA = California, USA. For acronyms of the location of each specimen see Table 1. 


jackknife test (by percentage contribution to the generation of 
the model) were annual temperature range {Pc — 58), 
precipitation in the driest month (Pc = 12.8) and mean diurnal 
range {Pc = 8.2). 

DISCUSSION 

Phylogeography and population genetic structure, -The 
observed mitochondrial genetic differentiation of Pardosa 
sierra populations from California versus the rest of the Baja 
California Peninsula coincides with the occurrence of the 
alleged vicariance event known as the northern-gulf (Riddle et 
al. 2000). Similar phylogenetic patterns have also been 
detected in this region in several vertebrate species, e.g., the 
toad Bufo punctatus (Riddle et al. 2000), the rodent Thomomys 
bottae (Alvarez-Castaheda & Patton 2004), the lizard Xantusia 
sp. (Sinclair et al. 2004) and the boa Lichanura trivirgata 


(Wood et al. 2008). The concordance in phylogeographic 
signal among these different taxa suggests that they may have 
responded in parallel to the same underlying process/es 
(Arbogast & Kenagy 2001; Zink 2002a). Therefore, it is 
probable that in this region, the northern-gulf vicariance event 
has influenced the genetic architecture of P. sierra. Both clades 
are reciprocally monophyletic (Avise 2000), and this is most 
likely due to the occurrence of an ancient barrier (Zink 2002a). 
However, low migration rate and reinforcement could be 
alternative explanations for this phylogeographic pattern 
(Munguia-Vega 2011; Dolby et al. 2015). Another explanation 
for reciprocal monophyly between clades would be that clade 
CA constitutes a cryptic species. The morphological charac¬ 
teristics of Californian specimens from this study (color 
pattern, size and male pedipalp) do not differ from those 
reported for P. sierra (Correa-Ramirez et al. 2010), however 
the female epigynum shows some variation in its internal 
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Figure 5.—Niche-based distribution models of Pardosa sierra, (a) Model for the last interglacial period (140,000-120,000 years ago), (b) 
Community Climate System Model corresponding to the last glacial maximum (LGM) 21,000 years ago. (c) Model for Interdisciplinary 
Research on Climate during the LGM. (d) Actual time model. Colored bars indicate the probability of a favorable habitat for P. sierra , with red 
indicating the highest value. Points in actual model (d) indicate the locations used to generate ENMs, and circles (b and c) indicate potential 
refugia during the LGM according to AUC = 0.75-0.99 range. 


anatomy (i.e., structure of the spermathecae and copulatory 
ducts); the latter provides some support for a possible cryptic 
complex, with diversification during the Late Miocene/Early 
Pliocene (Starrett & Hedin 2006) in P. sierra. Indeed, cryptic 
speciation in spiders has been detected in the same region 
(Ramirez & Chi 2004), and is potentially explained by high 
rates of diversification in spiders more generally (Starrett & 
Hedin 2006). To better delimit species we need to include more 
individuals and additional molecular markers (i.e., unlinked 
nuclear genes), perform detailed morphological analyses, and 
reconcile multiple gene trees (Starret & Hedin 2006) using new 
statistical methods based as coalescent-theory (Yang & 
Rannala 2010; Ence & Carstens 2011; Fujita et al. 2012). 

Although a phylogeographic break was observed around 
the mid-peninsular channel and Isthmus of La Paz (Crews & 
Hedin 2006) in spiders of the genus Homalonychus Marx, 1891 
(Homalonychidae), in our study phylogeographic breaks were 
not detected among P. sierra populations along the BCP (Fig. 
3). These results and a paucity of geological evidence refute the 
hypothesis of the occurrence of a mid-peninsular channel 
(Gristlier 2002; Murphy & Aguirre-Leon 2002; Crews & Hedin 
2006; Garrick et al. 2009), or if it did exist, it did not have a 
lasting impact on the genetic structure of P. sierra in this 


region. Indeed, although the signal from the COl data did not 
reveal population expansion, the wide distribution of a single 
haplotype (Hapl), the few haplotypes with restricted distri¬ 
bution (Fig. 2b) and the divergence time inferred (Fig. 4) all 
suggest that P. sierra populations were already present along 
the BCP at least 2.4 Mya. 

Fossil and rate calibrations for divergence date estimation 
using BEAST resulted in unsurprisingly different estimates of 
divergence times. We used a relaxed molecular clock model for 
both analyses, with the fossil calibrated analysis estimating 
much deeper ages for all nodes. Each method has its 
limitations, and using substitution rates alone has been 
criticized due to relaxation of the strict molecular clock in 
many taxa (Ho & Larson 2006; Ho 2007), and the problems 
posed by applying substitution rates from other species (Zink 
2002b). Meanwhile, the use of fossils to calibrate trees can be 
erroneous given the uncertainty of geological dating (Mag- 
allon 2004) and the incomplete nature of the geological record. 
In this study, we focus on the results of our fossil calibrated 
BEAST analysis, which are concordant with phylogeographic 
splits inferred for other vertebrate taxa (i.e., Peromyscus 
eren liens , Chaetodipns bailey /, Bufo pm wt at as) in the northern- 
gulf around 3 Mya (Riddle et al. 2000). 
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The genetic structure of P. sierra , as detected by microsat¬ 
ellite markers (Fig. 2c, d; Table 2), was not concordant with 
geographic breaks expected under a model of vicariance along 
the BCP. Despite the difference in mutation rates between 
microsatellite molecular markers and COl (Hedrick 1999), in 
this study both markers failed to detect phylogeographic 
signals related to the mid-peninsular and Isthmus of La Paz 
channels. This evidence could support the hypothesis that 
climatic phenomena configured the current genetic architec¬ 
ture of P. sierra populations, or that these phenomena erased 
any mitochondrial genetic signal that might have originated 
from these vicariance events. 

Pleistocene refugia.—The location of Pleistocene refugia 
depends on how each species responded to the cyclic climate 
oscillations of this period (Harl et al. 2014). This is evident in 
the few studies that have focused on identifying potential 
refugia during the Pleistocene in the Baja California Peninsula. 
For some plant species, such as totem pole cactus (Lophocer- 
eus schott //; Nason et al. 2002) and cardon (Pachycereus 
pringlei ; Gutierrez 2015), and some animal species, such as the 
gnatcatcher Polioptila califoniica (Zink et al. 2000), refugia 
have been observed towards the south of the peninsula. In 
contrast, for some arthropod species, the peninsula's central 
region acted as a refugium (Wilson & Pitts 2012; Graham et al. 

2014) . In the case of P. sierra , two geographic areas located in 
the north and south end of the peninsula most likely acted as 
separate refugia (Fig. 5b, c), and these areas coincide with the 
refugia suggested by Hafner & Riddle (1997). According to the 
niche-based distribution modeling analysis, P. sierra probably 
occupied a greater geographic range along the BCP during the 
last interglacial period (140,000-120,000 years ago; Fig. 5a). 
However, according to the simulation results, when the 
temperature decreased during the LGM (21.000 years ago), 
its potential habitat became fragmented. The lack of a 
relationship between latitude and microsatellite heterozygosity 
and between latitude and allelic diversity supports the presence 
of two refugia, because if only one had existed in the south, 
more genetic diversity would be expected in this region and 
there would be a gradient of decreasing genetic diversity 
towards the north (as has been suggested for other species in 
the BCP; Riddle & Hafner 2006; Valdivia 2014; Gutierrez 

2015) . Moreover, the existence of two populations of origin in 
opposing regions of the BCP could be the reason that a 
distance isolation structure was not detected. 

The presence of greater haplotype diversity in the popula¬ 
tions of the north and absence of exclusive haplotypes in the 
south (Fig. 2b), suggests that southern populations originated 
from those in the north before the refugia fragmented them. 
This is consistent with northern haplotypes of CA ( Hap7 and 
Hap3 ), which are more divergent relative to the other 
haplotypes (Fig. 2a). The wide distribution of only a single 
haplotype ( Hapl) across the peninsula suggests that this could 
be considered an ancestral haplotype (Avise 2000; Cuenca et 
al. 2003), derived from northern populations. The fact that 
Hapl is present on opposite sides of the BCP (Fig. 2b) 
supports the hypothesis that P. sierra took refuge in these 
locations on the peninsula. Likewise, the bottleneck detected 
using microsatellite loci in populations I and II (Fig. 2c) could 
explain the low COl diversity in the southern populations, 
given that this region coincides geographically with the area 


where COl transitions to a single haplotype (Hapl) from the 
north toward the south of the peninsula. However the window 
of time detected by BOTTLENECK software using microsat¬ 
ellite loci is only about ~25-250 generations in accordance 
with Cornuet & Luikart (1996) and Luikart et al. (1998), 
therefore it is necessary to analyze genetic markers with low 
mutation rates to test if a bottleneck could be related to the 
middle peninsular vicariant event. Finally, the intolerance of 
P. sierra to cold and its sensitivity to changes in rainfall or 
humidity (Van Dyke & Lowrie 1975) are constraits that limit 
these spiders to preferentially inhabit areas surrounding water 
bodies and oases (Punzo & Farmer 2006; Correa-Ramirez 
2010; Jimenez et al. 2015). Interestingly, variation in sea levels 
could be related to changes in the variables that most 
influenced the potential niche of the species during the 
Pleistocene (e.g., temperature and humidity), and may 
therefore be implicated in refugia occurring in more than 
one area of the BCP. Certainly, if cold intolerance forced P. 
sierra to occupy those disparate areas (Fig. 5b, c), as is likely, 
then it is plausible that postglacial recolonization processes 
were accountable for erasing the genetic signal that may have 
been produced by other vicariance events in the central and 
southern regions of the Baja California Peninsula. 

Based on these data, the northern-gulf channel appears to 
be the only vicariance event that left a lasting phylogeographic 
signal in P. sierra. This conclusion is supported by strong 
genetic differentiation between California (USA) populations 
and those from the Baja California Peninsula. The genetic 
differentiation patterns revealed by mitochondrial sequences 
and microsatellites failed to detect genetic breaks in regions 
where the mid-peninsular and Isthmus of La Paz canals were 
presumably formed. Two refugia were detected at both ends of 
the Baja California peninsula during the Pleistocene, specif¬ 
ically between the last interglacial period (140,000-120,000 
years ago) and the last glacial maximum period (2L000 years 
ago). 
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Appendix 1.—GenBank accesion numbers of mitochondrial (COI) Appendix 2 — Records of Pardos a sierra used to generate potential 

marker. The species listed were used in phylogenetic analysis for refuges for the Interglacial and during and after the Last Glacial 
Lycosidae. Maximum periods. All data were collected for this study in the states 

—— -——————-—-————— — of Baja California and Baja California Sur in Mexico and California 

Species _ Gen Bunk Accesion number in the United States of America. 


Alopecosa moriutii 

Arctosa littoral is 

AB564729 

JQ280374 

Location 

Lat 

Long 




Geolycosa xera 

DQ151816 

USA, CA, Sn. Juan Creek 

33.639 

-117.422 

Hogna helluo 

JQ280373 

USA, CA, Sn. Juan Creek 1 

33.607 

-117.444 

Lycosa tarantula 

KC550668 

USA, CA, Sn. Juan Creek, LakeLand Ville 

33.607 

-117.444 

Parclosa astrigera 

AY836055 

USA, CA, Margarita Truck Trail 

33.454 

-117.377 

Pardosa atromedia 

FJ 546466 

USA, CA, Sn. Felipe Creek 

33.066 

-116.553 

Pardosa atromedia 

FJ546467 

USA, CA, Descanso Town 

32.843 

-116.605 

Pardosa bellona 

JQ280362 

USA, CA, Pine Creek 

32.838 

-116.538 

Pardosa californica 

JQ280357 

Mexico, BC, Arroyo Las Palomas 

32.374 

-116.355 

Pardosa hamifeva 

JQ280364 

Mexico, BC, Arroyo Sn. Antonio Minas 

31.969 

-116.659 

Pardosa mi Ivina 

JQ280356 

Mexico. BC, Arroyo Sn. Salvado 

31.853 

-116.078 

Pardosa sierra Hapl 

FJ546464 

Mexico. BC, Arroyo Sn. Carlos I 

31.793 

116.501 

Pardosa sierra Hap2 

JQ280371 

Mexico, BC, Arroyo Sn. Carlos II 

31.786 

- 116.505 

Pardosa sierra Hap3 

JQ280366 

Mexico. BC, Ensenada 

31.783 

-116.600 

Pardosa sierra Hap4 

Fj546465 

Mexico. BC, Rancho Las Liebres 

31.584 

116.033 

Pardosa sierra Hap5 

KT364484 

Mexico, BC, Arroyo El Mejin 

30.980 

-116.095 

Pardosa sierra Hap6 

JQ280372 

Mexico, BC, El Rosarito 

28.617 

-114.033 

Pardosa sierra Hap7 

JQ280367 

Mexico, BCS, Sn Ignacio 

27.175 

-112.869 

Pardosa steva 

FJ546470 

Mexico. BCS. Cadeje 

26.367 

112.500 

Pardosa steva 

FJ546471 

Mexico. BCS, Carambuche-Sn. Isidro 

26.237 

112.002 

Pardosa sura 

FJ546468 

Mexico. BCS, Sn. lsidro-La Purfsima 

26.200 

-112.033 

Pardosa sura 

FJ546469 

Mexico. BCS, Arroyo Sn. Jose 

26.059 

-111.820 

Pardosa sura 

JQ280358 

Mexico, BCS. Arroyo Sn. Javier 

25.871 

111.546 

Pardosa vadosa 

FJ 546472 

Mexico. BCS, Sn. Pedro de la Presa 

24.833 

-110.983 

Pardosa vadosa 

FJ546473 

Mexico. BCS, El Pilar-Las Pocitas 

24.472 

-111.001 

Pardosa valens 

FJ546474 

Mexico, BCS. Rancho Camaron 

24.320 

-110.669 

Pardosa valens 

FJ 546475 

Mexico. BCS, Presa de la Buena Mujer 

24.088 

-110.191 

Pardosa sp 

JQ280361 

Mexico. BCS, Playitas 

23.986 

-110.187 

Pardosa sp 

JQ280360 

Mexico. BCS, El Novillo 

23.917 

-110.217 

Pardosa sp 

JQ280363 

Mexico, BCS, Camino a Sierra de La Laguna 

23.550 

-109.984 

Pirata canadensis 

KF368671 

Mexico, BCS, El Chorro, Santiago 

23.439 

-109.804 

Rabidosa rabida 

Trochosa ruricola 

DQ029232 
AB564731 

Mexico, BCS, Sierra de la Laguna 

23.233 

-109.867 




Venatrix pseudospeciosa 

JQ240195 












